Local Temperature and Universal Heat Conduction in FPU chains 
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It is shown numerically that for Fermi Pasta Ulam (FPU) chains with alternating masses and 
heat baths at slightly different temperatures at the ends, the local temperature (LT) on small scales 
behaves paradoxically in steady state. This expands the long established problem of equilibration 
of FPU chains. A well-behaved LT appears to be achieved for equal mass chains; the thermal 
conductivity is shown to diverge with chain length N as N 1 ^ 3 , relevant for the much debated 
question of the universality of one dimensional heat conduction. The reason why earlier simulations 
have obtained systematically higher exponents is explained. 



PACS numbers: 

It has long been established [l[ that Fermi Pasta Ulam 
(FPU) chains (one dimensional chains of particles with 
anharmonic forces between them) may not be able to 
achieve thermal equilibrium with seemingly reasonable 
initial conditions. This has interesting implications for 
the ergodic hypothesis, at the foundation of statistical 
mechanics. The results of Ref. fil led to the discovery of 
solitons in continuum versions [2[ , and eventually an un- 
derstanding of chaotic dynamics. With regard to the ini- 
tial results, a vast body of work has established [1,13 that 
quasi-periodic solutions exist below an energy threshold, 
above which the system does equilibrate. 

Although all this work has been with closed bound- 
ary conditions, the study of FPU chains has expanded 
to include heat bath boundary conditions: with slightly 
different temperatures imposed at the two ends, the ther- 
mal conductivity can be measured and shows anomalous 
properties [H, Q . Since the baths at the ends are at dif- 
ferent temperatures, the concept of equilibrium has to be 
extended: a local temperature (LT) that varies smoothly 
along the chain has to be defined. Surprisingly, this has 
not been fully investigated, even though the discussion of 
heat conductivity is in terms of Fourier's law [5[ which is 
meaningless if the local temperature is ill-behaved. When 
the heat baths at the two ends have equal temperatures, 
one can prove analytically that the onlypossible steady 
state is the one in thermal equilibrium [7[ ■ 

In this paper, we demonstrate for the first time that, 
for FPU chains with heat bath boundary conditions, the 
LT behaves paradoxically on small scales. In particular, 
we show through numerical simulations that for FPU- 
(3 chains with alternating light and heavy masses, con- 
nected to heat baths at temperatures Tr, and Tr at the 
left and right end respectively (with AT = Tr, — Tr suf- 
ficiently small that the system is in the linear response 
regime), the kinetic temperature of the particles oscil- 
lates as one moves along the chain. The ratio of ampli- 
tude of these oscillations to AT/N is constant as AT is 
reduced, and increases with the chain length N: in the 
vicinity of the 2i'th particle, the temperature difference 
between heavy and light particles scales approximately as 
ST ~ [AT / V~N]f (i/N) . Thus if one were to coarse grain 



over a region of the order of the mean free path, the intra- 
cell temperature uncertainty 0(AT/y~N) would domi- 
nate the 0(AT/N) change in temperature between adja- 
cent cells. (Similar results are obtained when the heavy 
and light particles are randomly ordered, so that no 0(1) 
coarse-graining length would work.) Nor is this a bound- 
ary effect with a characteristic decay length, since the 
dependence on i is through i /N. We also show that there 
are no problems when the FPU chain has equal masses. 
For one dimensional gases, lack of energy equipartition 
between heavy and light particles has been observed for 
hard particle systems^, but we have verified that T(x) 
as a function of position x (instead of particle number i) 
is smooth. For the FPU system, where the lattice con- 
stant can be taken to be arbitrarily large, this resolution 
of the problem is not applicable. 

For equal mass FPU chains, having verified the exis- 
tence of a well-behaved LT, we measure the heat conduc- 
tivity k(N) and demonstrate that 



k(N) ~ N a 



(1) 



with a = 1/3, in agreement with the earlier renormaliza- 
tion group (RG) analytical result [9( for a fluid model. 
We do this by simulating very long FPU chains with up 
to N — 65536 particles, showing that a = 1/3 is attained 
in this regime. This result is insensitive to system param- 
eters, in contrast to the apparent exponents for smaller 
N. We also explain why the apparent a for smaller N is 
systematically higher than 1/3, as seen in various earlier 
numerical simulations 0, [lO, [Tl|, dill]. This supports 
the assertion that there is only one universality class for 
heat conduction in one dimensional momentum conserv- 
ing systems (that have a well-behaved LT), contrary to 
earlier suggestions [l(| EH, E2, EH • 

The generalized FPU chain consists of a sequence of 
particles connected by springs between nearest neigh- 
bors. The Hamiltonian is J2i m i^i /2 + V( x i ~ x i+i)i 
where the potential energy of the interparticle springs 
is V(z) = k 2 z 2 /2 + k 3 z 3 /3 + k 4 z 4 /4 + ... . For all fig- 
ures shown in this paper, we use the FPU-/3 model where 
&2 = &4 = 1 and other fej's are zero. As noted in the orig- 
inal FPU paper [l[ , the incidental even symmetry of this 
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FIG. 2: (Color online) Kinetic temperature profile for a FPU- 
P chain with TV = 16384, T L = 2.0, T R = 0.5 and 7 = 2.0. 
The first three even moments of the velocity are shown; their 
agreement indicates a Gaussian velocity distribution. (Inset) 
Normalized temperature profiles for different TV. 



FIG. 1: (Color online) Oscillations of the kinetic temperature 
profile for FPU-/3 dimer chains with mass ratio 2.62, 7 = 
2. (Top inset) Full profile for a system with TV = 128 and 
T L = 2.0, T R = 0.5. (Top) TV = 128 particles, T L , R = 2 ± 
AT/2. The scaled temperature is (T - 2)/AT; the different 
plots line up. (Bottom) (T 2i - T 2l +i)N 1/2 versus 2i/N for 
different TV. Tl,r = 2.0,0.5. The TV 1/2 and 2i/N are chosen 
to approximately match the vertical and horizontal scales of 
the plots. Together, the figures imply ST « [AT/VN]f(i/N). 



potential may cause non-mixing of even and odd modes, 
but this is not the case with the heat bath boundary con- 
ditions used here. After fixing AT, we allow the system 
to reach steady state and then measure the kinetic tem- 
perature Ti = mi(vf) of each particle. Unless otherwise 
noted, the end particles i = 1 and i — TV are connected to 
Langevin baths at temperatures Tl — 2.0 and Tr = 0.5 
respectively, by adding damping and noise terms to their 
equations of motion, which are integrated with an accu- 
rate Verlet-like algorithm p^ |. 

FigQ] shows Ti as a function of particle number i for a 
dimer chain with TV = 128 and mass ratio 2.62. An un- 
usual oscillating temperature profile is seen; the lighter 
particles are hotter than the heavier ones on the left and 
colder on the right. Since the sign, not just the magni- 
tude, of 2j_|_i — Ti oscillates, it is not because of a thermal 
conductivity k that oscillates on the microscopic scale. 
The behavior shown here is robust to changes in the mass 
ratio, although the oscillation amplitude changes, or the 
interparticle potential (e.g. if an exponential potential is 
used, i.e. a Toda lattice with alternating masses). Os- 
cillations in Ti as a function of i are also seen when the 
Langevin baths are replaced with Nose Hoover baths [l5| . 

Figfl] also shows that if AT is reduced or TV is in- 
creased, the temperature difference between the 2i'th 
and (2i + l)'th particles scales approximately as ST ~ 



[AT/y/N]f(i/N). Thus coarse graining over a region of 
size O(l), comparable to the mean free path, creates an 
unusual LT: the uncertainty in temperature in a coarse- 
grained region is greater than the variation between ad- 
jacent regions. Moreover, the 'decay length' over which 
the oscillations penetrate into the interior of the chain 
is a fixed fraction of TV, showing that this is not negligi- 
ble for large TV. Non-monotonic temperature profiles are 
also seen when the heavy and light particles are ordered 
randomly, in which case the heavy and light kinetic tem- 
peratures track two separate smooth curves. 

Our results show that anharmonicity and disorder are 
not sufficient for a well-behaved LT even with heat bath 
boundaries imposing a 0(1/TV) temperature gradient. 
When a well-behaved LT is achieved, e.g. (as we will 
show) for equal mass FPU chains, it should be viewed as 
being fragile. The question remains: what are the neces- 
sary and sufficient conditions for a local temperature? 

We also simulate FPU-/3 chains of equal mass particles 
(of unit mass) as described in the previous paragraphs. 
When steady state is reached, we observe an approxi- 
mately linear temperature profile as shown by Fig. 
There is a slight curvature near the boundaries, which 
decreases with TV. Fig.© also shows that the velocity 
distributions are Gaussian (at least for large TV) which 
is necessary for a LT. This indicates that the equal mass 
FPU chain with heat baths at the boundaries has a well- 
behaved LT. 

With this reassurance, we proceed to measure the heat 
conductivity as a function of TV for the equal mass case. 
With a small temperature difference applied across the 
system, Fourier's law predicts that the heat current j 
should be equal to — kVT, with a K that depends on 
microscopic properties. Equivalently, j — — «( AT/TV), 
where we have defined VT as AT/TV. By measuring 
j(N), deviations from j ~ 1/TV are interpreted as a TV- 
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dependent conductivity n(N) and a consequent break- 
down of Fourier's law. 

For one dimensional momentum conserving systems 
where a well-behaved LT exists, a RG study of the hy- 
drodynamic equations of a normal fluid d showed that 
Eq.(TTJ is satisfied with a = 1/3. This has been confirmed 
by simulations of hard particle gases fill [T7L El , al- 
though very large systems are required [191 ] and the issue 
is not completely settled [13] • On the other hand, nu- 
merical simulations of oscillator chains, including FPU 
chains, give various exponents @, El, El, El, El for 
different systems, often slightly higher than 1/3. This 
seems consistent with early results from mode-coupling 
theory (MCT), whi ch p redict a heat conductivity ex- 
ponent of a = 2/5 [la EE El, although recent MCT 
analyses predict exponents that depend on the leading 
nonlinearity [H, HlJ and the extent of transverse mo- 
tion El . The apparent agreement between the numerical 
and MCT results has led to speculation that there may 
be two (or more) universality classes with different expo- 
nents [ll El- The most recent MCT analysis [H, |2l| 
predicts that a ^ 1/3 is restricted to even potentials 
V(z), which is why we have studied the FPU-/? model 
here. 

Recently, two of us have extended the earlier RG treat- 
ment to systems with broken symmetry [23 ] . which can 
occur on intermediate length scales in one dimension (as 
in nanotubes). The resulting crystalline hydrodynamic 
equations also yielded a = 1/3, like the earlier fluid re- 
sult [11. From numerical results, it was argued that 
the apparent a > 1/3 found for FPU chains is probably 
a crossover effect from hard particle systems. However, 
the possibility that the crossover could be pushed out to 
N — ► 00 could not be ruled out. Thus the numerical 
and analytical evidence so far allow for FPU chains to 
be a singular limit for the RG with a special value of 
the conductivity exponent a, motivating our numerical 
simulations. 

The heat current flowing in steady state is measured 
as a function of N for equal mass FPU chains. The time 
averaged current j is defined by j = — XiV'(xi+i — 
Xi)/N), where Xi is the displacement from equilibrium of 
the i th particle. As shown in Fig. [3J n(N) = -jN/AT 
satisfies Eq.fll) for large N with a = 0.333 ± 0.004, in 
strong agreement with the RG prediction. 

To test the sensitivity of this result to different baths, 
we run simulations with Langevin baths with different 
damping constants 7=0.4, 2, and 10. We also replace 
the stochastic Langevin baths with deterministic Nose- 
Hoover [151 ] thermostats, for which we use the fourth or- 
der Runge-Kutta integrator. Fig.(j3j) compares the RG 
prediction and the MCT prediction for systems with 
these different baths and bath parameters. As can be 
seen in the figure, an asymptotic exponent of 1/3 is at- 
tained for all these systems, whereas the apparent ex- 
ponents for smaller N depend on system parameters. 
Moreover, it is possible to understand the deviation of 
the apparent exponent from 1/3 for small system sizes. 
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FIG. 3: (Color online) Heat current as a function of N. 
Root mean square errors from O(10 5 — 10 6 ) measurements are 
shown except when they are smaller than the points. (Top) 
Conductivity versus N. The last five points fit to a slope of 
a = 0.333 ± 0.004. The baths are described in Fig.©. (Bot- 
tom) jN 1-01 versus N for a = 1/3 and a = 2/5. In the large 
N regime, a is definitely less than 2/5 and appears to agree 
quite well with the 1/3 prediction for all data sets. Langevin 
baths with 7=0.4, 2, and 10, and one data set with Nose- 
Hoover baths, are shown. 



As shown in Ref.[6(, if the damping constant for the 
Langevin baths is very large or small, there is a large 
'contact resistance' at the boundaries of the chain. The 
current only depends weakly on N, resulting in an ap- 
parent a > 1/ 3. ( Similar considerations apply to Nose- 
Hoover baths [2J].) This is confirmed by our results: 
the plot for 7 = 2 reaches the asymptotic limit fastest, 
whereas 7 = 0.4, 10 have apparent exponents closer to 
0.4 for small N. Taken together, the large-iV exponent 
of 1/3, the universality with respect to bath parameters, 
and the explanation for how the apparent exponent be- 
haves as a function of N and 7 for small N convincingly 
supports the RG prediction. 

Since Fourier's law is only applicable in the linear re- 
sponse regime, we halve the temperature difference be- 
tween the ends and simulate the same spring system with 
Langevin baths with 7 = 2, T L = 1.625 and T R = 0.875. 
Though the error bars are larger, a still agrees with 1/3, 
verifying that the system is in the linear response regime. 

The exponent a measured in our simulations of FPU 
chains clearly differs from the measurements from other 
simulations [H, El, El, El • This disagreement is mainly 
because very large system sizes are needed. Moreover, 
we use a step size h — 0.0025 — 0.005 that is an order of 
magnitude smaller than the step size used for the dynam- 
ics in (ll]; by comparing numerical and exact results for 
harmonic springs, we have found that a small h is neces- 
sary for proper convergence. Finally, we have shown the 
results as a function of the coupling to the bath, which 
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has not been done before, and explained why the baths 
increase the apparent a, especially if the coupling is not 
optimal. We note that only equilibrium correlations for 
a periodic FPU-/3 chain are shown in Ref. PI]. As dis- 
cussed in the second paper in Ref. [Hj], this is problem- 
atic. 

However, in Ref. pl| . where purely quartic springs 
(/c2 = 0) and system sizes similar to our simulations 
were used, a ~ 0.44 was extrapolated from an indi- 
rect measurement and a > 0.4 was measured directly 
from nonequilibrium simulations. We have found a — 
0.38 — 0.39 for this system. Because the analytical re- 
sults all obtain the same asymptotic a for the FPU-/3 and 
quartic models, and we have shown that the apparent a 
decreases with N for the FPU-/3 model, we believe that 
a for the quartic system will also asymptotically change 
to 1/3. The alternative, that a for the FPU-/3 chain will 
reverse its change with TV and revert to ~ 0.4, seems un- 
likely. Indeed, the latest MCT results [2l[ obtain a = 0.5 
for even potentials (including the quartic model and the 
FPU-/3 model) which is quite far from earlier [l3| and 
our numerical results. Nevertheless, it is unclear why the 
pure quartic system should need exceptionally large N. A 



final resolution of the issue requires an analytical demon- 
stration of an error in one of the competing methods. 

The existence of a paradoxical LT on small scales for 
slight alterations to the equal mass FPU chain could be 
partly responsible for the unusually large N needed to 
reach the asymptotic limit (most notably for purely quar- 
tic potentials). Conversely, slight changes can improve 
the conv erg ence to a = 1/3 as seen by adding a cubic 
term [H,[ll(, transverse motion [l2j], or collisions fli"ll22l ]. 
Recent simulations of nanotubes obtain a = 1/3 [251 ]. 

In summary, we have shown that the local tempera- 
ture behaves paradoxically when an 0{\/N) temperature 
gradient is applied to FPU chains with unequal masses; 
coarse-graining with any O(l) averaging length does not 
cure this. This renders existing analytical models for heat 
conductivity inapplicable. However, a well-behaved LT 
is established for equal mass chains, or when the baths 
are at the same temperature Q. For equal mass chains, 
large scale simulations support the much debated heat 
conductivity exponent being 1/3 as predicted 0. 

We thank Sriram Ramaswamy, Ram Ramaswamy, and 
Peter Young for useful comments. 
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